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ABSTRACT 

We present a refined gravitational lens model of the four-image lens system B1608-I-656 based on new 
and improved observational constraints: (i) the three independent time-delays and flux-ratios from Very 
Large Array (VLA) observations, (ii) the radio-image positions from Very Large Baseline Array (VLBA) 
observations, (iii) the shape of the deconvolved Einstein Ring from optical and infrared Hubble Space 
Telescope (HST) images, (iv) the extinction-corrected lens-galaxy centroids and structural parameters, 
and (v) a stellar velocity dispersion, dap = 247 ± 35 km s~^, of the primary lens galaxy (Gl), obtained 
from an echelle spectrum taken with the Keck-II telescope. The lens mass model consists of two elliptical 
mass distributions with power-law density profiles and an external shear, totaling 22 free parameters, 
including the density slopes which are the key parameters to determine the value of Hq from lens time 
delays. This has required the development of a new lens code that is highly optimized for speed. The 
minimum-x^ model reproduces all observations very well, including the stellar velocity dispersion and 
the shape of the Einstein Ring. A combined gravitational-lens and stellar dynamical analysis leads to a 
value of the Hubble Constant of Ho=75l^ km s'^Mpc"^ (68% CL; rJm=0.3, Oa=0.7). The non-linear 
error analysis includes correlations between all free parameters, in particular the density slopes of Gl and 
G2, yielding an accurate determination of the random error on Hq. The lens galaxy Gl is '^5 times more 
massive than the secondary lens galaxy (G2), and has a mass density slope of Jq-^ = 2.03lg :[4±0.03 (68% 

CL) for p (X , very close to isothermal {j'—2). After extinction-correction, Gl exhibits a smooth 
surface brightness distribution with an i?^/'* profile and no apparent evidence for tidal disruption by 
interactions with G2. Given the scope of the observational constraints and the gravitational-lens models, 
as well as the careful corrections to the data, we believe this value of Ho to be little affected by known 
systematic errors (^5%). 

Subject headings: gravitational lensing — cosmology: distance scale — galaxies: elliptical and 
lenticular, cD — galaxies: structure — galaxies: kinematics and dynamics 



1. INTRODUCTION 

The physics behind gravitational lensing is well under- 
stood and primarily based on gravity, i.e.. General Rela- 
tivity (e.g., Schneider, Ehlers, & Falco 1992). Since Refs- 
dal (1964) it has been known that arcsecond-scale strong 
gravitational lenses - those with multiply-imaged sources 
- provide a tool to measure the expansion rate of the Uni- 
verse, i.e., the Hubble Constant (Ho), if the mass distri- 
bution of the lens(es) and the time-delay(s) between the 
lensed images are known. This provides an elegant one- 
step and global method to measure Ho, independent of 
local distance-ladder techniques (e.g., Parodi et al. 2000; 
Saha et al. 2001; Freedman et al. 2001), which could be 
prone to the accumulation of unknown systematics, and 
seem difficult to improve beyond the current 10% preci- 
sion on the value of Hq . 

More recently, CMB observations with WMAP have 
been combined with different cosmological data sets (i.e., 
Ly-a and Large Scale Structure observations) to indi- 
rectly infer a value of iio—7lt.^ km s^^Mpc^^ (Bennett 



et al. 2003; Spergel et al. 2003) similar to that from 
the Hubble Space Telescope (HST) Key-Project, who find 
Ho=72+8 km s'^Mpc"^ (Freedman et al. 2001). Because 
Ho can not be measured directly from the CMB data alone, 
a combination with other datasets is necessary, and thus it 
is prone to its own systematic errors (e.g. Seljak, McDon- 
ald, & Makarov 2003; Bridle, Lahav, Ostriker, & Stein- 
hardt 2003). Even so, the agreement between the dif- 
ferent values of Hq from different techniques, including 
those from Sunyaev-Zeldovich observations (e.g.. Mason, 
Myers, & Readhead 2001), is striking. Since Hq is one 
of the most fundamental cosmological parameters, how- 
ever, its measurement should be subject to multiple cross- 
examinations, with the ultimate goal of breaking the 10% 
precision limit. 

Although gravitational lensing can provide a straightfor- 
ward and potentially very precise measurement of Ho, its 
major problems are that it requires (i) accurately mea- 
sured time-delays obtained from long monitoring cam- 
paigns and (ii) an accurate determination of the mass dis- 
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Fig. 1.— A compilation of the HST-F160W images of B1608+656: (a) The original reduced image of B1608+656. (b) The extinction- 
corrected image (see text), (c) The model surface brightness distributions of Gl and G2, fitted to the extinction-corrected brightness 
distributions, (d) The original image of B1608+656 after subtraction of the lens galaxy models extincted by the inferred dust-screen. 



tribution in the lenses and the field (i.e., the "lens poten- 
tial"). At present about a dozen lens systems have mea- 
sured time-delays with difi^erent degrees of precision (see 
Courbin, Saha. & Schcchtcr 2002, for a recent summary). 
This has shifted attention to understanding the often un- 
derestimated complexity of the mass distribution of the 
lens. For example, the first discovered gravitational lens 
system Q0957+561 (Walsh et al. 1979), has an exquisitely 
measured time-delay with an error ^1% (e.g. Kundic et al. 
1997; Oscoz et al. 1997; Haarsma, Hewitt, Lehar, & Burke 
1999; Colley & Schild 2000; Oscoz et al. 2001; Slavcheva- 
Mihova, Oknyanskij, & Mihov 2001; Goicoechea 2002; 
Ovaldsen, Teuber, Schild, & Stabell 2003), but a complex 
lens potential that includes a nearby cluster. Even though 
much efi^ort has been put into constraining this particular 
lens system, the lens potential is still not known to ade- 
quate precision to allow a satisfactory determination of Hq 
(Keeton et al. 2000, and references therein). 
One of the dominant degeneracies in gravitational lens 



models - even for systems that are relatively isolated and 
have no major external perturbers (e.g., groups or clusters) 
- is that between the radial mass profile and the inferred 
value of Hq. Often point-image constraints (i.e., positions 
and fiux-ratios) can be satisfied equally well by mass dis- 
tributions that have very different slopes of their radial 
mass profile (e.g., Wambsganss & Paczynski 1994; Refs- 
dal & Surdej 1994; Witt, Mao, & Schechter 1995; Witt, 
Mao, & Keeton 2000; Wucknitz & Refsdal 2001; Surpi & 
Blandford 2001; Wucknitz 2002; Williams & Saha 2000). 
This is similar to the well-known mass-sheet degeneracy 
(Gorenstein, Shapiro, & Falco 1988). In general, steeper 
(shallower) mass profiles lead to higher (lower) inferred 
values of Hq . To break this degeneracy and reliably mea- 
sure Ho, one has to determine the mass profile of the lens 
galaxies. Often the most effective way is to use additional 
external constraints, for example from stellar kinematics of 
the lens galaxy (e.g., Treu & Koopmans 2002a; Koopmans 
& Treu 2002). 
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An often cited result is that of the gravitational 
lens system PGl 115+080, which gives a value of 
Ho=44±4 kms-^Mpc"^ (Impey et al. 1998) for an 
isothermal lens mass distribution {p oc r~^). Not only is 
this value low compared to most other lens and non-lens 
measurements, but the quoted error does not include one 
of the major sources of uncertainty - the unknown radial 
mass profile. In fact, early-type galaxies can have a consid- 
erable scatter in their radial mass profiles (e.g., Gerhard et 
al. 2001) around the effective radius, and steeper mass den- 
sity profiles can increase Hq well above 60 km s~^ Mpc~^ 
(Impey et al. 1998). 

It was recently shown that, if one includes the ob- 
served stellar velocity dispersion (Tonry 1998) of the lens 
galaxy in PG1115+080 as additional constraint, the ef- 
fective slope of the mass density profile can be measured 
directly and also properly included in the error estimate on 
Hq. The effective slope of the density profile for this partic- 
ular system is found to be steeper than isothermal, leading 
to a higher value of Ho=59t7^ ± 3 km s"^ Mpc"^ (Treu & 
Koopmans 2002b). This measurement is in better agree- 
ment with most other methods and its larger error bars 
include the contribution of the residual uncertainty on the 
mass distribution of the lens. This analysis has led us to 
initiate a new program with Keck and the Very Large Tele- 
scope (VLT) to measure the stellar velocity dispersion of a 
number of additional systems with known time-delays, in 
order to determine a more precise value of Hq from grav- 
itational lensing, less affected by the radial mass profile 
degeneracy. 

In this paper, we focus on the gravitational lens sys- 
tem B1608+656 (source redshift Zg — 1.39; lens redshift 
zi = 0.63) (Myers et al. 1995; Snellen et al. 1995; Fassnacht 
et al. 1996). The system is unique in that all three time de- 
lays are known between its four lensed images (Fassnacht 
et al. 1999, 2002b), with accuracies of a few percent. In ad- 
dition, multi-color images are available in the HST archive 
that show the lensed arcs of the radio-source host galaxy 
connecting to a full Einstein Ring (e.g., Blandford, Surpi, 
& Kundic 2001; Kochanek et al. 2001; Surpi & Bland- 
ford 2003). The multi-color images allow us to correct the 
positions and surface brightness distributions of the lens 
galaxies for extinction, removing a major source of system- 
atic error common to previous modeling efforts. Finally, 
we have measured the stellar velocity dispersion of the 
dominant lens galaxy from a spectrum obtained with the 
Echelle Spectrograph and Imager (ESI; Sheinis et al. 2002) 
on the Keck-II Telescope. With this additional informa- 
tion and the improved constraints on the time-delays, we 
are well-suited to refine the lens models of this system and 
remove many of the major degeneracies and systematic 
biases in the determination of Hq. 

This paper is organized as follows. In §2 and §3 we dis- 
cuss the observations and data reduction of the HST and 
Keck data, respectively. In §4, we describe the observa- 
tional constraints obtained from VLA and VLBA radio 
observations'^, HST images, and the Keck spectrum. In 
§5, we briefly discuss a new lensing code and the min- 

^ The National Radio Astronomy Observatory is a facility of the 
Associated Universities, Inc. 

^ IRAF (Image Reduction and Analysis Facility) is distributed by t 
Association of Universities for Research in Astronomy under cooper; 



imization procedure. In §6 and §7, the modeling results 
and inferred value of Hq are discussed. In §8, we summa- 
rize and discuss our results. Throughout this paper, we 
assume fim = 0.3 and = 0.7. 

2. IMAGING 

The main objectives of the analysis of the multi-color 

HST imaging are to obtain: (i) accurate extinction- 
corrected images of the two main lens galaxies, from which 
their true centroids and structural parameters (i.e., posi- 
tion angle, ellipticity and effective radiiis) can be deter- 
mined; and (ii) a deconvolved Einstein Ring after sub- 
traction of a model of the lens galaxies. The shape of 
the Einstein Ring places additional constraints on the az- 
imuthal structure of the lens model (Kochanek et al. 2001). 
We note that the analysis of the HST images presented 
here is independent of that in Surpi & Blandford (2003; 
hereafter SB03), but produces nearly identical color maps. 
Finally, we bring all astrometric data obtained from the 
HST images to a common "reference frame", i.e., that of 
the high-resolution VLBA observations. This allows us 
to combine the radio and optical/infrared observations in 
an self-consistent way, taking proper care of potential dif- 
ferences (i.e., offsets, rotations and scalings) between the 
different data sets. 

2.1. Observations and data reduction 

The HST images of the gravitational lens system 
B1608+656 are available from the HST archive. Since 
the primary concern here is obtaining high spatial reso- 
lution and an accurate model of the point-spread function 
(PSF), we select the images with the best resolution and 
sampling in each of the available bands (see also the dis- 
cussion in SB03): 41,344s in 7 exposures with the Near 
Infrared Camera and Multi Object Spectrograph (NIC- 
MOS) Camera 1 (NICl) through filter F160W (GO-7422; 
PI: Readhead); and 4 exposures on the Planetary Camera 
(PC) of the Wide Field and Planetary Camera 2 (WFPC2) 
through two filters for total exposure times of 11,500s 
and 11,600s in the F814W and F606W bands, respectively 
(HST GO-6555; PI: Schechter). 

The images are first reduced using a series of iraf^ 
scripts based on the drizzle package (Fruchter & Hook 
2002). The scripts clean and combine the images, by it- 
eratively refining the measurement of relative offsets and 
the identification of cosmic-rays and defects (similar to the 
scripts used in Koopmans & Treu 2002, Treu et al. 2003). 
The final F160W "drizzled" image is shown in the top left 
panel of Fig.l. Images of the system in the other filters are 
similar to the ones shown in SB03 and are not repeated 
here. 

2.2. Dust extinction correction 

The B1608-I-656 system is characterized by two lens 
galaxies, Gl and G2, with a strong color gradient in the 
region connecting the two galaxies. The positions of the 
centroids of the two galaxies with respect to the multi- 
ple images were observed to vary by as much as 86 milli- 

National Science Foundation operated under cooperative agreement by 

le National Optical Astronomy Observatories, which are operated by the 
i,tive agreement with the National Science Foundation. 
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arcseconds (mas) between the F160W and F606W im- 
ages (SB03; Koopmans & Fassnacht 1999; Fassnacht et 
al. 2002b). 

To remove this source of systematic uncertainty, we as- 
sume that the observed color gradients arise from dust 
extinction (see SB03) and that the two lens galaxies have 
the same intrinsic colors. We use the available colors to 
measure the color excesses and correct for dust extinction 
based on an extinction law (see below for a discussion). 
To check for systcmatics, we perform the correction for 
each pair of filters, in the reference frame of each one of 
the filters. The results from each pair of filters are very 
similar and affect the results only at a very small level, 
which will be taken into account in the final error analy- 
sis on Hq. In the interest of space we will describe as an 
example only the procedure for the F814W-F160W pair 
in the F160W reference frame, mentioning the differences 
in the results for the other frames where relevant. First, 
we deconvolve the F814W image using a Lucy Richardson 
algorithm and a PSF generated using tiny tim (Krist 
1993). Second, we align the deconvolved F814W image 
to the F160W image assuming that the positions of the 
point sources are the same in each filter and not affected 
by dust'^, and that the transformation is a rotation and 
translation with arbitrary axis scaling. Third, the decon- 
volved and transformed F814W image is convolved with 
a NICMOS PSF generated by TINY TIM, to match the 
resolution of the F160W image. Fourth, the ratio of the 
matched images in the two filters is used to produce a 
color map, which is found to be very similar to the one 
shown in SB03. Fifth, it is assumed that the intrinsic col- 
ors of Gl and G2 are spatially uniform and equal to the 
minimum of the color measured in the region including 
Gl and G2 (see also Surpi & Blandford 2003). Finally, 
the color excess is converted into extinction in the indi- 
vidual filters using the following relations computed using 
the Galactic extinction law by Cardelli et al. (1989) for 
Rv = 3.1: A(F160W)/Av=0.415, A(F814W)/Av=1.091, 
A(F606W)/Av=1.561. The extinction-corrected image^ 
in F160W is shown in the top right panel of Fig. 1. 

This procedure is repeated for each pair of filters, in 
the reference frame of each filter. We emphasize that this 
procedure does not assume a particular smooth surface 
brightness profile for the lens galaxies, nor does it force 
the centroids of the lenses to coincide in the different fil- 
ters. However, the production of a smooth extinction- 
corrected surface brightness distribution for the lens galax- 
ies strongly suggests, although does not guarantee, that 
the extinction correction has been done properly. Simi- 
larly, the consistency of the astrometry between various 
filters provides an additional test of the uncertainties re- 
lated to the assumption of a specific dust extinction law, 
the assumption of coincident centroids of the lensed im- 
ages, and the PSF modeling used for convolution and de- 
convolution purposes. These tests are discussed in the 
next two sub-sections, together with surface photometry 
and astrometry measurements. 



Surface photometry is performed on the extinction- 
corrected F814W and F160W images, as described in 
Treu et al. (1999, 2001). The F606W image is not used 
for this task, given that it has the largest dust extinc- 
tion and the lowest signal-to-noise ratio (SNR), especially 
on G2. The galaxies are modeled as two R^^^ profiles 
- which has been shown to be a good description of 
the extinction-corrected surface brightness profile of Gl 
(Blandford, Surpi, & Kundic 2001) - that are fit simul- 
taneously to obtain effective radii (Re) and total magni- 
tudes (Table 1). Only the region well inside the Einstein 
Ring is used in the fit, to minimize contamination by the 
extended ring structure. Very similar photometric param- 
eters are obtained in the two filters, demonstrating that 
uncertainties in the alignment of the images and in the 
PSF modeling do not affect significantly the derived sur- 
face photometry. The best fitting model in filter F160W 
is shown in the lower left panel of Fig. 1 . 

To inspect the residuals from the fit, we apply the dust 
extinction map to the best fitting i?^/^ model and subtract 
the reddened model from the original image. The bright- 
ness contribution of the galaxies at the ring radius is found 
to be negligible. The residual image is shown in the lower 
right panel of Fig.l, demonstrating that indeed the data 
are consistent with two relatively smooth spheroids, red- 
dened by a dust lane (Surpi & Blandford 2003) without 
strong evidence that the primary lens galaxy Gl is much 
affected by tidal interactions with G2. The ring structure 
is very clear in the residual image (e.g., Kochanek et al. 
2001; Surpi & Blandford 2003), and is even more promi- 
nent in the deconvolved residual image (Fig. 5). 



Gl 



G2 



F160W (mag) 
F8UW (mag) 
Re,Fi60W (arcsec) 
Re,F8i4w (arcsec) 
b/a={l - e) 
Major axis P.A.(°) 



16.41±0.08 
18.23±0.12 
0.62±0.05 
0.58±0.06 
0.5±0.1 
79 ±2 



18.18±0.05 
19.99±0.10 
0.27±0.03 
0.26±0.04 

= 1 

N/A 



Table 1 — Observed photometric quantities of the lens galaxies 
Gl and G2. All magnitudes are in the Vega system. The sys- 
tematic uncertainty on the NICMOS zero point (see the NIC- 
MOS instrument handbook; Mobashcr 2002) is not included. 
Magnitudes arc tiic total magnitudes obtained from fitting the 
galaxy surface- brightness profiles with an i?^''* law. Magni- 
tudes have been corrected for internal extinction, as discussed 
in the text, but are not corrected for Galactic extinction. The 
axial ratio of G2 is set to unity during the fit. 



2.4. Astrometry 

The extinction-corrected images arc used to measure the 
locations (in pixels) of the centroids of Gl and G2, using 
the task imcentroid in iraf. The deconvolved residual 
images are used to measure the centroids of the multiple 
images of the lensed source, as well as the shape of the 
ring (see §4.2 and right panels in Fig. 5). 

Since we perform all measurements from the HST im- 
ages in their original pixel frames, we need to align the 
optical and infrared frames to the much more accurate 

^ Since the centroids of the lensed images are very sharp, they change very httle due to extinction (see also §2.4). 
* Note that this correction is only made to the lens galaxies and not to the Einstein Ring. 



2.3. Surface photometry 
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VLBA radio frame. We determine the transformation (i.e., 
rotation, translation and scaling) that minimizes the posi- 
tional differences between the lensed image centroids mea- 
sured from the HST images, and their radio counterparts, 
weighting properly by the measurement errors. We allow 
for different pixel scales in the x/y directions on the NICl 
and PC chips. We find pixel scales of 0'.'0457/0'.'0456 for 
the PC and 0'.'0434/0'.'0431 for NICl, in excellent agree- 
ment with the values in the HST instrument handbook 
(Mobashcr 2002). 

The positional differences between the optical/infrared 
and radio image positions are < 10 mas in RA and DEC 
for each of the individual lensed images. Since the trans- 
formation is determined from the full set of four images, 
we estimate an error of ~4 mas in x and y on the trans- 
formation between the radio and NICMOS frames. Using 
this transformation, we map the trace of the brightness 
peak of the Einstein Ring to the radio frame (see §4.2). 
Since the transformation errors are negligible compared 
with the positional determination of the brightness peak 
along spokes, one can assume for all practical purposes 
that the radio data and the Einstein Rings are perfectly 
aligned. 

We use the same transformation to bring the galaxy cen- 
troids in F160W to the radio frame. The combined cen- 
troid errors and transformation errors quadratically add 
to an error of 6 mas in x and y. We adopt this value as 
the centroiding errors of the galaxies with respect to the 
lensed radio images. The measurement errors on the radio 
positions are negligible {<^1 mas; Koopmans & Fassnacht 
1999). We repeat the same procedure for the F814W im- 
age. 

After transformation to the radio frame, the positions of 
the lens galaxies Gl and G2 agree to within 14 mas and 20 
mas between the F160W and F814W filters, respectively, 
as opposed to 72 mas and 38 mas in the uncorrected frame. 
Slightly larger differences are found using the F606W filter, 
as expected because this is the filter where uncertainties 
in the dust corrections are the most severe^. All astromet- 
ric measurements from the F160W and F814W filters are 
listed in Table 2. The deconvolved Einstein Rings in the 
F160W and F814W filters are shown in the left panels of 
Fig.5. 



3. SPECTROSCOPY 

In this section we describe the spectroscopic observa- 
tions of B1608-I-656 with ESI. The main objective of the 

observations was to obtain the stellar velocity dispersion of 
the primary lens galaxy (Gl), the measurement of which 
allows us to determine the power-law slope of its radial 
mass profile and break the degeneracy with the inferred 
value of the Hubble Constant. Special care is required 
in order to determine an unbiased stellar velocity disper- 
sion and a correct uncertainty, since Gl is a post-starburst 
galaxy. This is accomplished through extensive Monte 
Carlo simulations of artificial spectra with different frac- 
tions of "old" and "young" stellar populations. 

^ The maximum difTerence in the centroid coordinates between F606W 
to 86 and 45 mas respectively in the uncorrected images of SB03. 

^ developed by D. Sand and T. Treu; Sand et al. (2003), in prep. 
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Fig. 2. — Spectrum of B1608-I-656— Gl in one of the cchelle or- 
ders around the G— band (in counts, not flux calibrated). The lower 
histogram shows the noise per pixel. A boxcar smoothed version of 
the spectrum is shown above the original spectrum to guide the eye. 
The wavelength range used for the kinematic fit was 6750-7200A 
(observed), which includes the G-band, H7, and Fe4384 features, 
but excludes H(5. The atmospheric B-band (hatched band) and H7 
were masked out in the kinematic fit. 



3.1. Observations and data reduction 

We observed B1608-h656 using ESI (Sheinis et al. 2002) 
on the Keck-II Telescope on 2000 July 3, for a total in- 
tegration time of 5400s (3xl800s). The seeing was good 
with FWHM« 0'.'85 and the night was clear. Between 
each exposure, we dithered along the slit to allow for a 
better removal of sky residuals at the red end of the spec- 
trum. The slit (20" in length) was positioned at PA= 83°, 
i.e., aligned to within 4 degrees of the major axis of Gl. 
The slit- width of 0'.'75 yields an instrumental resolution 
of fT~20km s~^ which is adequate for measuring the stel- 
lar velocity dispersion and removing narrow sky emission 
lines. The centering of the galaxy in the slit was con- 
stantly monitored by means of the ESI viewing camera - 
the galaxy was bright enough to be visible in exposures 
of a few seconds duration - and we estimate the centering 
perpendicular to the slit to be accurate to better than O'.'l. 

Data reduction was performed using the iraf package 
EASI2DS as discussed in Koopmans & Treu (2002) and 
Sand et al. (2002,3). A one dimensional spectrum was ex- 
tracted (Fig. 2) by summing the signal within an aperture 
corresponding to 1'.'7 x 0'.'75. The SNR of the spectrum in 
the vicinity of the G-band (~4304A) is 20 A" ^ Note the 
very prominent Balmer absorption features H7 and H(5, 
typical of a K+A post starburst population (Dressier & 
Gunn 1983; Myers et al. 1995; Surpi & Blandford 2003). 

and F160W is now 21 mas, with an rms scatter of 14 mas, as opposed 
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Table 2 

astrometry of the b 1608+656 system 



VLBA/8.4-GHZ HST/F814W HST/F160W 



Image 


Ax (") 


Ay (") 


Ax (") 


Ay (") 


Ax (") 


Ay (") 


A 


=0.0000 


=0.0000 


-0.0003 


+0.0009 


-0.0029 


+0.0042 


B 


-0.7380 


-1.9612 


-0.7376 


-1.9588 


-0.7342 


-1.9498 


C 


-0.7446 


-0.4537 


-0.7455 


-0.4528 


-0.7469 


-0.4446 


D 


+1.1284 


-1.2565 


+1.1298 


-1.2563 


+1.1365 


-1.2546 


Gl 


N/A 


N/A 


+0.4124 


-1.0592 


+0.4261 


-1.0581 


02 


N/A 


N/A 


-0.3091 


-0.9306 


-0.2897 


-0.9243 



3.2. The Stellar velocity dispersion of Gl 

The luminosity weighted velocity dispersion Cap of Gl 
was measured with the Gauss-Hermite Pixel-Fitting Soft- 
ware (van der Marel 1994) on the spectral region covering 
the observed wavelength ~ 6750 — 7200 (Fig. 2) excluding 
B.S and including the G-band, H7, and Fe4383 (Trager et 
al. 1998). As kinematic templates we used spectra of G K 
giants observed at twilight with a 0'.'3 slit width, appro- 
priately smoothed to match the instrumental resolution 
of the 0'.'75 slit. The regions around 6870 A, affected by 
atmospheric absorption, and the H7 feature were masked 
during the fit. The result is robust with respect to changes 
in the spectral region used in the fit, the fit parameters 
(e.g., the order of polynomial used to reproduce the con- 
tinuum level), and the adopted template. The best fit 
gives: a^p ~ 247 ± 35 km s~-^ in an aperture of l'.'7x0'.'75 
centered on Gl. 

To minimize the amount of corrections to the data, we 
will only use this measurement of the luminosity weighted 
velocity dispersion (as opposed to the central velocity dis- 
persion, defined in a standard aperture) in the kinematic 
analysis (§4), properly taking the seeing and aperture ef- 
fects into account by applying those corrections to the 
models. However, to facilitate comparison with earlier 
work and other applications of this measurement, we also 
give the central velocity dispersion^, a = 269 + 40 km s"""^, 
within a standard circular aperture of radius Rc/8. 

We note here that if one adopts the measured central ve- 
locity dispersion and the structural parameters measured 
from the extinction corrected images (Table 1), to deter- 
mine the offset of Gl from the local Fundamental Plane, 
one finds the galaxy to be significantly brighter than quies- 
cent early- type of the same mass and at similar redshifts. 
This is consistent with what is typically found for K+A 
galaxies (e.g., van Dokkum & Stanford 2003). In contrast, 
this result differs from the findings of Rusin et al. (2003a), 
who estimated a larger velocity dispersion from the im- 
age separation and measured a lower luminosity for Gl, 
resulting in an offset from the local Fundamental Plane 
consistent with quiescent evolution. 

The rest of this section describes Monte Carlo simula- 
tions aimed at ensuring that our measurement provides 
not only an unbiased value of the stellar velocity disper- 
sion of Gl, but also a realistic estimate of the uncertainty. 



In fact, although the SNR and resolution are more than 
adequate to measure the velocity dispersion of a uniformly 
old stellar population (e.g., Treu et al. 2001), the spectrum 
of Gl is more challenging, because of the contamination 
by young stars. Our Monte Carlo approach works as fol- 
lows. First, we model the spectrum as a combination of a 
G8III star, to represent the old population, and of an A5V 
star, to represent the young population (Dressier & Ounn 
1983; Franx 1993). The stellar spectra are smoothed to 
a velocity dispersion of 250 km s^^, taking into account 
the initial spectral resolution of the stellar spectra*, and 
then combined with variable weights. Second, for each 
set of weights (hereafter "old" fraction and "young" frac- 
tion) 1,000 realizations arc constructed, adding noise to 
reproduce the SNR of our data. Third, we run the Gauss- 
Hermite Pixel-Fitting Software (van der Marel 1994) on 
each realization to simulate the stellar velocity dispersion 
measurement. 



;r~-3oo 



100 



0.2 0.4 0,6 
young fraction 



0.4 0,6 0,8 1 

line strength 



Fig. 3. — Monte Carlo simulations of the stellar velocity disper- 
sion recovery from artificial spectral of B 1608+656— Gl. The dotted 
horizontal line shows the input value, the pentagons shows the av- 
erage value measured from 1,000 realizations. The heavy error bars 
show the average uncertainty estimated by the code, while the light 
error bars show the r.m.s. scatter of the recovered values 

The results of the Monte Carlo simulations are shown 
in Fig. 3. In the left panel we plot the average and r.m.s. 
scatter (solid pentagons with thin error bars) of the recov- 
ered velocity dispersion rrnt as a ftmction of the fraction 
of A5 V light (young fraction). For young fractions up to 
~ 0.6, the procedure recovers the input velocity dispersion 
to within 4%. Only for a young fraction larger than 0.7 
does the procedure become biased. This happens because, 
as the young fraction increases, the information content 



^ The correction from luminosity weighted velocity dispersion to central velocity dispersion is computed based on the range of observed velocity 
dispersion profiles of local E/SO galaxies, and taking into account seeing effects as described in Treu et al. (2001). 

* The G8III star is from the ESI data, the A5 V is taken from the library of Jacoby, Hunter & Christen 1984; see Treu et al. 2001 for a 
discussion of the resolution of that library. 
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from the metal absorption features gets washed away by 
the contribution of the young population. Similarly, for 
young fractions less than 0.6, the r.m.s. scatter (thin error 
bars) agrees remarkably well with the average uncertainty 
returned by the code (thick error bars). We can estimate 
the fraction of young stars in the spectrum of Gl by com- 
paring the observed depth of the G-band and of H7 with 
those in the model spectra. In Fig. 4 we show that a com- 
posite spectrum with ~40% of young population matches 
the observed spectrum, while fractions ^70% and higher 
produce much deeper absorption lines than observed. 
We therefore conclude that the fraction of young stellar 
light in Gl must be smaller than ^70% and that the value 
of Gap is the correct and unbiased measurement of the stel- 
lar velocity dispersion of Gl. Similarly the error estimate 
is unbiased. 



4200 4300 4400 




6800 6900 7000 7100 7200 



wavelength (A) 

Fig. 4.— A fit of o-ap for Gl in B1608+656 for two difl^'erent frac- 
tions of the young stellar population (shifted for clarity). A young 
fraction of 70% clearly overpredicts the depth of the H7 feature. 

A more quantitative confirmation of this result comes 
from the inspection of the line strength, which is used by 
the Gauss-Hermite Pixel-Fitting Software (see definition 
in van der Marel 1994) to measure the depth of the absorp- 
tion features: the larger the line strength, the deeper the 
absorption features. In the models, by construction, the 
line strength of the old stellar population anticorrelates 
very well with the young fraction. In fact a very similar 
conclusion is drawn if line strength is used as independent 
variable to show the result of the Monte Carlo simula- 
tions (right panel of Fig. 3). When the line strength of 
the old stellar population is larger than 0.4 (correspond- 
ing approximately to a young fraction smaller than 0.7) 
the signal to noise ratio and resolution of our spectrum 
are sufficient to obtain an unbiased measurement of the 
stellar velocity dispersion. For a real composite popula- 
tion the interpretation of this parameter is more complex, 
since it depends on the chemical abundances and possible 

^ Some early-type galaxies show evidence for iso-surface-density twists 



additional contributions to the continuum. Nevertheless, 
it provides a useful quantitative measure of the depth of 
the metal absorption features, which is the main parame- 
ter that controls the accuracy of the fit. The best fitting 
template yields line strength 0.5 ± 0.05, also confirming 
that the measurement of velocity dispersion is unbiased. 

4. MASS MODEL AND OBSERVATIONAL CONSTRAINTS 

We model the Icnsing and kinematic properties of both 
lens galaxies, Gl and G2, assuming that their luminous 
plus dark-matter mass distributions are elliptical (i.e., 
oblate in three dimensions'' ) and have a radial power-law 
dependence, (i.e., p oc r^'' ). In particular, we assume for 
the lens models a dimensionless surface mass density of 

n{x,y) = h[x' + {y/qeff-''^'\ (1) 

where is the axial ratio of the surface mass distribution 

and 7' the slope of the corresponding density distribution 
(7' = 2 is designated as "isothermal" henceforth). The 
mass distribution has its centroid at {xe,ye) and a major- 
axis position angle of 9^ (measured north through east). 
In addition, we allow for a single external shear for the 
system with strength 7cxt and position angle 0cxt • 

We choose to use power-law mass models, since they 
have been successful in reproducing many of the de- 
tailed lensing observations (e.g., Kochanek 1995; Cohn, 
Kochanek, McLeod, & Keeton 2001; Munoz, Kochanek, & 
Keeton 2001; Rusin et al. 2002; Winn, Rusin, & Kochanek 
2003) and more recently also kinematic observations (e.g., 
Treu & Koopmans 2002a; Koopmans & Treu 2003) of 
lens systems and lens ensembles (e.g., Rusin et al. 2003a), 
even though we acknowledge that this choice might not 
be unique. We feel at this point, however, that there is 
no strong incentive to go beyond these models. From a 
practical standpoint, including the constraints from the 
Einstein Ring becomes computationally very expensive for 
more complex lens models (see §5 for a discussion). 

For the dynamical modeling (i.e., determining the stel- 
lar velocity dispersion), we solve the spherical Jeans equa- 
tions, assuming a spherical mass distribution with the 
same power-law index 7' as used for the lens models (i.e., 
Eq.(l)). We model the luminous mass distribution as a 
trace population embedded in the luminous plus dark- 
matter potential - with a Hernquist luminosity- density 
profile (Hernquist 1990). We have also examined Jaffe 
(1983) models and find differences typically <1% for the 
models that we explore. Henceforth, we only use the Hern- 
quist profile, which closely follows a R^/'^ brightness pro- 
file over the radial range of interest and is analytically 
tractable. Since the break radius of the dark-matter halo 
is well beyond the Einstein radius, we assume it to be 
much larger than the latter, with negligible effects on the 
calculation of kinematic quantities (see Treu & Koopmans 
2002a; Koopmans & Treu 2003, for a full discussion). Fi- 
nally, we allow for anisotropy of the stellar velocity ellip- 
soid, modeling the radial anisotropy parameter /3 (Binney 
& Tremaine 1987) with a Osipkov-Merritt (Osipkov 1979; 
Merritt 1985a, b) model with anisotropy radius ri. For 
Vi = OG, the stellar velocity dispersion is fully isotropic. 
We also take the seeing and aperture into account. We 

that could be due to triaxial mass distributions. 
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refer to Koopmans & Treu (2003) for a further discussion 
of our kinematic model. 

We now continue with a more detailed description of 
the precise set of constraints that we use in the lensing 
and kinematic models that are presented in §6. 

4.1. The Radio and Optical Image Constraints 

The positions of the four lensed radio images can most 
accurately been determined from high-resolution radio ob- 
servations. We use the VLB A positions listed in Koop- 
mans & Fassnacht (1999), which are consistent with VLA 
observations to within ~1 mas (Fassnacht ct al. 2002b). 
We choose positional errors on the VLBA image; positions 
of 1 mas, which are considerably larger than the formal 
measurement errors (tens of micro-arcseconds). This er- 
ror is approximately that by which the most massive mass 
substructures in the lens (expected to be on the order of 
~ 10^ Mq) could potentially change the image positions. 
Since mass substructures are "randomly" distributed such 
positional shifts can practically be regarded as random er- 
rors. 

Radio flux ratios might not be reliable constraints on 
macro lens models; either because of mass substructure 
(e.g. Mao & Schneider 1998; Mctcalf & Madau 2001; 
Keeton 2001; Chiba 2002; Metcalf & Zhao 2002; Dalai 
& Kochanek 2002; Bradac et al. 2002; Keeton et al. 2003; 
Bradac et al. 2003), radio microlensing (e.g., Koopmans & 
de Bruyn 2000; Schechter & Wambsganss 2002) or inter- 
stellar-medium (ISM) propagation effects (Koopmans et 
al. 2003). Hence, we conservatively adopt a large error 
of 20% on the individual radio fluxes (i.e., -^28% on the 
flux ratios). This is supported by the fact that the formal 
errors on the flux ratios within a single season are much 
smaller than the differences between seasons (Fassnacht ct 
al. 2002b). For the flux ratios we adopt the approximate 
median value from three independent monitoring seasons 
of B1608-I-656 (Fassnacht et al. 2002b). Since previous 
models had difficulties reproducing the flux ratio of image 
D (Koopmans & Fassnacht 1999; Fassnacht et al. 2002b), 
it is not used as a constraint. Hence, if the refined lens 
models can recover the flux of image D, we would have 
additional confidence in the lens models. 

Finally, we use the time delays between the VLA images 
as constraints (Fassnacht et al. 1999, 2002b). Time-delay 
ratios are particularly strong constraints on the mass pro- 
file, since they are almost a direct measure of the depth of 
the potential at the lensed-image positions. Since there are 
three independent time delays and one Hubble Constant 
to solve for, they add two constraints on the mass model 
(we emphasize here that time delays are affected by sub- 
structure only at a negligible level, <Cl%). All point-image 
constraints are listed in Table 3. 

4.2. The Einstein-Ring Constraints 

The F160W image of B1608+656 exhibits a fully- 
connected Einstein Ring around the two lens galaxies 
(Kochanek et al. 2001; Surpi & Blandford 2003). The 
structure of Einstein Rings can be used as additional con- 
straint on the gravitational lens models (Blandford, Surpi, 
& Kundic 2001; Kochanek et al. 2001). In particular, the 
unit vector along a radial spoke at the position of the 
brightness peak (on the spoke), when projected on the 



source plane, is perpendicular to the gradient of the surface 
brightness distribution of the unlensed source (Kochanek 
et al. 2001). When the isophotes of the source are as- 
sumed to be elliptical - often a reasonable assumption - 
the trace of brightness peaks (along radial spokes) places 
extra constraints on the lens potential at the minimal cost 
of adding four additional free parameters (ring-source po- 
sition, ellipticity and position angle). Moreover, for an 
elliptical source, the trace of the ring is independent of as- 
smnptions on the radial behavior of the surface brightness 
distribution of the source. 

Even though the SNR decreases with increasing distance 
from the central nucleus of the lensed source along the ring, 
it is still sufficient to fully trace the ring (see Fig.l and 
Kochanek et al. 2001). This can be done even more accu- 
rately on the deconvolved image, where the Einstein Ring 
is "sharper" in the radial direction (Fig. 5). We note that 
color gradients over the image have very little effect on 
the brightness peaks of the sharp ring (see also Kochanek 
et al. 2001). To trace the deconvolved Einstein Ring, we 
define 90 independent spokes, separated by 4 degrees in 
angle (consistent with the improved image resolution after 
deconvolution) , radiating from a point roughly in the mid- 
dle of the ring. The origin of the spokes is arbitrary. The 
brightness distribution of the deconvolved ring along each 
spoke is subsequently fitted with a Gaussian. From this 
fit, the position of the brightness peak and the FWHM of 
the ring are obtained. The resulting fitted ring is shown 
in Fig. 5. To ensure that the ring is not significantly af- 
fected by uncertainties produced by the deconvolution, we 
repeat the same procedure for the F814W image. We as- 
sume that the error on the ring peak is proportional to 
the ring width, since it is hard to estimate an error from a 
deconvolved image. We adopt 0.42 xctp as error for F160W 
(with fTr=FWHM/2.35 from the Gaussian fit), which gives 
X^/DOF^l, only weakly dependent on the particular lens 
model. This particular choice is not crucial, however, and 
somewhat higher or lower scale factors can be chosen with- 
out affecting the modeling results. Hence, we put more 
weight on those parts of the ring that are sharpest (near 
the radio images) and less weight on the less well-defined 
parts of the ring. Notice the remarkable agreement be- 
tween the traces of the ring in F160W and F814W (Fig.5), 
lending credit to the trace (compare also to Kochanek et 
al. 2001, who find a very similar trace from the original, 
not deconvolved, ring). 

4.3. Constraints from the Lens Galaxies 

As for the Einstein Ring, we use the transformation dis- 
cussed in §2.4 to bring the positions of the two lens galax- 
ies in the F160W and F814W images to that of the VLBA 
frame with positional errors of 6 mas with respect to the 
radio images (Table 3). The good agreement of the galaxy 
centroids from different filters - after extinction correc- 
tion - seem to indicate that the extinction corrections ap- 
plied to the F160W image have been done properly. Con- 
fidence in the extinction correction is further supported by 
the fact that, after correction, Gl is well described by a 
i?^/^ brightness profile (see also Blandford et al. 2001). It 
seems highly unlikely that an improper extinction correc- 
tion would produce a well-behaved smooth surface bright- 
ness distribution from an originally patchy image. In the 
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Table 3 

Galaxy centroid positions and point-image constraints on the B1608-I-656 system 



Imag 




Ay (") 


*5'norm 


At (days) 


A 


=0.0000±0.001 


=0.0000±0.001 


2.020±0.404 


31.5±1.5 


B 


-0.7380±0.001 


-1.9612±0.001 


1.000±0.200 


=0.0 


C 


-0.7446±0.001 


-0.4537±0.001 


1.034±0.207 


36.0±1.5 


D 


+1.1284±0.001 


-1.2565±0.001 


0.347±oo 


77.0±1.5 


Gl 


+0.4261±0.006 


-1.0581±0.006 


N/A 


N/A 


G2 


-0.2897±0.006 


-0.9243±0.006 


N/A 


N/A 



"The positions for components A-D are from Koopmans & Fassnacht (1999). 
''The time delays and flux ratios are from Fassnacht et al. (2002b). 
°The Gl and G2 positions are from this paper. 



lens models, we therefore adopt the light centroids from 

the extinction-corrected F160W images, transformed to 
the radio frame, as the final values for the mass centers of 
Gl and G2. 

In addition, modeling of the extinction-corrected bright- 
ness distribution gives a position angle of (79±2)° (north 

to east) of Gl. Since the extinction-corrected galaxies ap- 
pear relatively smooth and unperturbed it is safe to as- 
sume that Gl is most likely not much affected by tidal 
interactions with G2. The position angle of Gl is used as 
a prior on the position angle of its mass distribution. 

As a final model constraint, we use the stellar veloc- 
ity dispersion of Gl, cxap = 247 ± 35 km s~^ (§3,2), which 
provides information on the slope of its radial mass profile. 

5. LENSING CODE AND MINIMIZATION 

To model the B1608-I-656 system, we have developed 
a new lens code, highly optimized for speed and accu- 
racy. This is necessary because most of the models require 
X^-minimization in a 22-dimensional parameter space (see 
§6.1 for a listing of all free model parameters). Since the 
structure of the Einstein Ring is included as a constraint, 
each iteration (i.e., calculation) is computationally ex- 
pensive, especially if there are no analytic solutions for the 
lens potential. Fiu'thermorc, a full non-linear error analy- 
sis of the parameters of the lens model requires a repetitive 
minimization of the function for different sets of fixed 
lens parameters. 

To maximize speed, the code is first optimized in its 
use of memory, which can be critical when copying large 
data structures (i.e., those containing all constraints and 
intermediate results). Second, analytic solutions are used 
wherever possible, as well as publicly-available fast algo- 
rithms developed for specific power-law mass models (e.g., 
Barkana 1998; Chae, Khersonsky, & Turnshek 1998; Chae 
2002). A significant speed-up of the x^ minimization and 
error analysis is gained by using the MINUIT minimization 
package from CERN (James 1998), which has been opti- 
mized for finding function minima in large-dimensional pa- 
rameter spaces with a minimum number of function evalu- 
ations. We have compared the results from the new code, 
for B1608-I-656 and various simple test systems, with those 
from a previously used code (e.g., Koopmans & Fassnacht 

The non-linear error analysis consists of re- minimizing for a range of fixed values for the parameter of interest (James 1994) . The 68% 
confidence limits, quoted in this paper, correspond to the parameter range with Ax^ = 1 w.r.t. the minimum-x^ solution. 



1999) and from the publicly-available code from Keeton 
(2001). We find that they all agree to within the expected 
numerical errors. 

The increase in speed of the new code over the other 
codes, when we include the Einstein Ring data for 
B1608-I-656, can be several orders of magnitude, de- 
pending on the specific situation. For example, x^^ 
minimization with 22 free parameters usually takes min- 
utes, compared to hours for the code by Keeton (2001), 
when we use mass models that are not analytically 
tractable (e.g., non-isothermal power-law density profiles). 
A full non-linear error analysis for all parameters a built- 
in option with MINUIT takes several hours with the new 
code, but only several minutes extra for a single parame- 
ter (e.g., Ho). A similar analysis with the code by Keeton 
(2001) is expected to take several days and is clearly not 
feasible for systems of this complexity. 

To find the global x^ minimum, we first optimize all 
free model parameters using only the radio data as con- 
straints. This rapidly leads to a model that fits all the 
radio constraints extremely well. In the second step, the 
Einstein Ring data are added to the constraints. Since the 
mass model parameters are already relatively close to the 
best solution, it takes only several minutes to converge to 
the x^ minimum. To check that indeed the correct global 
X^ minimum has been found, the models are started from 
different initial conditions and also using different mini- 
mization techniques (i.e., simulated annealing and down- 
hill simplex methods). Monte Carlo minimization is done 
once a x^ minimum is reached. Finally, the minimum in 
all parameter-space directions is traced while performing 
a non-linear error analysis^°. None of this has led us to 
suspect that the global minimum has not been found. 
This is also supported by the fact that the observed lens 
properties are reproduced well by the lens models. 



6. RESULTS 

In this section, we describe the results from the lens 
models. After listing the free parameters and adopted pri- 
ors of the general model in §6.1, we consider in §6.2 a 
simplified problem in which the slopes of the galaxies are 
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Table 4 

Lens mass model parameters for B1608+656. 
SIE SPLEl SPLE2 SPLEl+D SPLEl+D SPLE2+U SPLE2+D 











(iso) 


(aniso) 


(iso) 


(aniso) 


xe (arcsec) 


0.425 


0.426 


0.425 


0.425 


0.425 


0.425 


0.426 




-0.291 


-0.290 


-0.290 


-0.291 


-0.290 


-0.291 


-0.290 


yg (arcsec) 


— 1.U71 


-1.069 


— i.oby 


— i.oby 


— i.oby 


— i.oby 


— i.oby 


-0.929 


-0.928 


-0.928 


-0.928 


-0.928 


-0.928 


-0.928 


b (arcsec''' 


0.531 


0.553 


0.553 


0.526 


0.555 


0.531 


0.553 


0.288 


0.263 


0.263 


0.269 


0.263 


0.268 


0.263 


Qe 


0.606 


0.607 


0.606 


0.604 


0.603 


0.605 


0.606 




0.340 


0.308 


0.308 


0.318 


0.307 


0.316 


0.308 


( ) 


YD. 8 


77.0 


77.0 


77.0 


TT A 

77.0 


TT A 

77.0 


TT A 

77.0 


69.2 


68.5 


68.5 


68.4 


68.4 


68.4 


68.5 


7 


=z.uu 
=2.00 


1 QQ+0.20 
2.12 


1 00+014 
^•^^-0.15 

2.12 


9 nf^+014 

^■"°-0.13 

2.12 


9 nn+014 

^■^^-0.15 

2.12 


9 04+0.11 

2.12 


9 AA+0.12 

2.12 


7cxt 


0.085 


0.081 


0.081 


0.077 


0.081 


0.078 


0.081 


^'ext (°) 


6.8 


10.6 


10.6 


13.4 


10.4 


12.9 


10.6 


Xg (arcsec) 


0.095 


0.103 


0.102 


0.088 


0.102 


0.090 


0.102 


Us (arcsec) 


-1.070 


-1.070 


-1.070 


-1.069 


-1.070 


-1.069 


-1.070 


Qh 


0.639 


0.640 


0.640 


0.634 


0.640 


0.635 


0.640 


Oh n 


113.2 


113.0 


113.0 


112.1 


113.1 


112.3 


113.0 


Ho (km/s/Mpc) 


99.8 


74±i^ 
98.0 


74i^ 
98.0 


76±^ 
98.1 


74i^ 
98.0 


98.1 


98.0 



Note: The first (second) line for each pajrameter denotes values for Gl (G2). The 68% confidence limits on Hq and have been determined 
from a full non-lineaj: error analysis, where all free parameters are varied for a range of fixed values of these parameters until increases by 
unity. We limit ourself to these two most important parameters, since this calculation is computationally very expensive. 



fixed to isothermal. This is done to facilitate comparison 
with previous work and to illustrate the effect of correcting 
the centroid positions for dust extinction. In Section 6.3 
we explore the full range of parameter space, and hence its 
effects on Hq, allowing the slopes to vary to reproduce the 
lensing geometry. Finally, in Section 6,4 we add the dy- 
namical constraints to reduce the uncertainties and break 
the mass profile degeneracy of Gl. All results have been 
calculated using the new lens code and arc based only 
on constraints from the radio data, the NICMOS F160W 
images, and the stellar kinematics. Constraints from the 
WFPC2 F814W images will be used as additional check 
on systematic effects in Section 7.3. 

6.1. Free Parameters and Priors 

To allow maximum freedom in the lens models and a 

proper calculation of the significance of the derived value 
of the Hubble Constant, we allow nearly all of the model 
parameters (i.e., position, lens strength, position angle, el- 
lipticity, and density slope) of each lens galaxy to vary, 
constraining several with Gaussian priors. The source 
adds three additional free parameters (i.e., position and 
flux), the Hubble Constant adds one, and the external 
shear adds two. The model of the elliptical source that 
describes the Einstein Ring (§4.2) adds an additional four 
free parameters (position, ellipticity and position angle). 
In total there are up to 22 free parameters. 

We use Gaussian priors on the galaxy positions, placing 
them at their observed positions and allowing for 1-cr er- 



rors of 6 mas in x and y (Table 3). Hence, increases by 
unity if the galaxy positions deviate by l-tr from the ob- 
served positions. Similarly, we allow the position angle of 
Gl to vary using a Gaussian prior of 79±2 degrees. Since 
B1608+656 is part of a small group (e.g., Fassnacht et al. 
2002a) we allow for an external shear. The group is not 
massive and has a low velocity dispersion (Fassnacht et 
al. in prep.). We therefore place a conservative O.OOiO.lO 
(1-cr positive) prior on the strength of the external shear, 
although this prior is not particularly important. 

Finally we place a Gaussian prior on the density slope of 
the secondary galaxy (G2), for which we have no kinematic 
information. We choose 7^3 = 2.00 ± 0.10, fully consis- 
tent with the spread in values of 7' found from 0047-285 
and MG2016-F211 (Treu & Koopmans 2002a; Koopmans 
& Treu 2003). Since G2 is '^5 times less massive than Gl, 
its contribution to the properties of the lens system is less 
pronounced. We discuss this prior further in §7.3. 

6.2. Isothermal lens model 

To allow a comparison with previous modeling efforts 
(Koopmans & Fassnacht 1999; Kochanek et al. 2001; Fass- 
nacht et al. 2002b), we first set the density slopes of both 
Gl and G2 to isothermal (i.e. 7'=2). The structure of 
the ring is included as a constraint and the priors are as 
discussed in §6.1. 

We have listed the model parameters of the minimum- 
solution in Table 4, where this model is designated 
as SIE. Compared with previous isothermal lens models 




Fig. 5. — The deconvolved Einstein Rings of B1608+656 in the F160W (upper left) and F814W (lower left) bands, after extinction correction 
and subtraction of a model of the two lens galaxies. The contour levels increase by factors of two (arbitrary units). The solid curve indicates 
the outer critical curve (see text) of our best lens model (§6.4). The right panels indicate the trace of the Einstein Ring (thick solid line) and 
the ±1— 0"r width of the ring (two thin solid lines). To smooth the trace, we plot it for 360 spokes (only 90 are used in the modeling; see text). 
The dot-dashed curves show the minimum-x^ SPLEl+D models. 



(Koopmans & Fassnacht 1999; Fassnacht et al. 2002b), 
we find that the inferred value of the Hubble Constant, 
Ho=7i;^ kms-^Mpc-^ (68% C.L), has increased by 

about 6-8 km s'^ Mpc"\ The error includes all random 
errors and correlations between free parameters. This in- 
crease is partly due to the two main lens galaxies being 
closer to each other after extinction correction. Although 
each galaxy is isothermal, their combined potential be- 
comes steeper, leading to an increase in Hq (e.g., Wucknitz 
2002; Williams & Saha 2000). Other contributions come 



from including external shear and widening the errors on 
the image flux ratios and positions, which change the 
surface. 

Besides reproducing all point-image constraints well, an 
excellent fit (nearly identical to Fig. 5) to the boxy struc- 
ture of the Einstein Ring is also found. A comparison by 
eye to the fit in Kochanek et al. (2001; their Fig. 5) shows 
that our fit is significantly better, in particular between the 
image pairs B-C and D-A. The difference can probably be 
attributed to the fact that we use a deconvolved Einstein 
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Table 5 

Recovered point-image properties from the minimum-x^ SPLEl+D (isotropic) model. 



Image 


Ax (") 


Ay (") 


'S'norm 


At (days) 


A 


+0.00005 


+0.00003 


2.0408 


32.42 


B 


-0.73823 


-1.96125 


0.9313 


0.00 


C 


-0.74461 


-0.45379 


1.0869 


36.20 


D 


+1.12845 


-1.25630 


0.4230 


76.51 



Ring shape and extinction-corrected lens galaxy ccntroids. 
We find a source axial ratio of = 0.64+0.03 and position 
angle of Og = (113 ± 2)°. The axial ratio is similar, but 
the position angle of the ring source is significantly differ- 
ent from the value in Kochanek et al. (2001). The optical 
source has a minor offset from the source for the radio im- 
ages (<0'.'03) which could, for example, be attributed to a 
small asymmetry (e.g. lopsidedness or m=l mode) of the 
host galaxy or deviations from ellipticity. Forcing the cen- 
troid of the host galaxy to coincide with the radio source 
increases the of the ring fit (by ~30%), but changes the 
resulting mass model parameters only marginally. Since 
m=l modes for the brightness distribution of galaxies are 
not uncommon and automatically accounted for in models 
with a free host-galaxy centroid, we choose not to force 
such a coincidence. We find values identical to those de- 
rived from the F160W image when the ring shape from the 
F814W data is used, and also when using the code from 
Keeton (2001). 

As a further check, we calculate the stellar velocity dis- 
persion of Gl based on this isothermal lens model. De- 
tails of the calculation method are outlined in Koopmans 
& Treu (2003) and §4. We also come back to this in more 
detail in §6.4. We find (Tap — 230 km s~^ for isotropic mod- 
els (rj = oo) and trap = 250 km for anisotropic models 
with n = i?eff, inside the proper aperture and corrected 
for seeing effects. These values agree remarkably well with 
the measurement of dap = 247 + 35 km s^^, even though 
this information has not been used in the isothermal mod- 
els. This suggests that the simplest isothermal lens model 
is probably not too far from the final model. 

6.3. Power-law lens models 

The next step is to allow the slope of the dominant lens 
galaxy (Gl) to vary. Since the isothermal models show 
that Gl is ^5 times more massive than G2, the contribu- 
tion of G2 to the structure of the Einstein Ring and the 
lensed images is relatively small. In general we find that 
the properties of G2 are ill-defined if we leave its density 
slope (or ellipticity; see §7.3) entirely unconstrained. We 
therefore place a prior on its slope, = 2.00 ± 0.10, 
as suggested from similar mass models of other lens sys- 
tems (§6.1). The slope of Gl is left fully unconstrained. 
We designate this model as the "Singular Power-Law El- 
lipsoid" (SPLE) mass model. The SPLEl model is that 
model with a prior only on the power-law slope of G2. 

The results are listed in Table 4. We find a slightly bet- 
ter fit to the data (i.e. Ax^=— 1-8) than for the isother- 
mal models. Most notably, the slope of Gl is found 
to be = I.99I0J0' nearly identical to isothermal. 

We note that the effect of changing the aperture, effective radius, 



The slope of G2 is somewhat steeper than isothermal, 
which leads to a slight increase in the Hubble Constant, 
Ho=74tJ? km s'^ Mpc"^ (68% C.L) compared with the 
SIE mass models. Since the density slopes are included as 
free parameters, the resulting error on Hq is larger than 
for the isothermal models which have fixed density slopes. 
The error on Hq is less than 15%, even without kinematic 
constraints. If we include a similar prior on the slope of Gl 
as on G2 (the SPLE2 model), the results arc again nearly 
identical (Table 4). The error on Hq decreases consider- 
ably, because it is dominated by the uncertainty on the 
density slopes, as discussed previously. As for the isother- 
mal models, the point-image constraints and the shape of 
the Einstein Ring are reproduced well. 

6.4. Lensing & Dynamical models 

The final step in the modeling effort is to include the 

stellar velocity dispersion of Gl as an additional con- 
straint. This is not as simple as for a case with a single 
lens galaxy. Since we have two lens galaxies, we can not 
simply attribute all of the mass within the Einstein radius 
(or critical curve) to only the lens galaxy Gl. 

We therefore follow a different approach. First we deter- 
mine the minimum-x^ parameters for a range of models 
with fixed values of 7qi ranging from 1.6 to 2.4. We then 
calculate the mass enclosed within an aperture with a ra- 
dius of 1" centered on Gl, using the lens model parame- 
ters. Using the enclosed mass and the slope of the density 
profile of Gl, we can then determine the stellar velocity 
dispersion inside the observed aperture (for details, see §4 
and Koopmans & Treu 2003). The model velocity disper- 
sion is compared to the observed value and the resulting 
Ax^) assuming a Gaussian error, is added to the x| from 
the lens constraints alone. 

Since h, q and 7' all enter into the determination of the 
mass enclosed within 1", the minimum of xf + Ax^ is not 
necessarily the true minimum. However, since Yq-^ is by 
far the parameter that most dominates the determination 
of the inferred stellar velocity dispersion^ ^, for all prac- 
tical purposes the difference between this approach and 
the correct approach (which is to calculate the dispersion 
at every x^ evaluation) is negligible compared to the final 
error on Hq. We note that calculating the dispersion at ev- 
ery minimization step would slow down the minimization 
process by several orders of magnitude and is not feasible. 

The resulting best-fit models for the lensing and dy- 
namical models, designated as SPLEl+D and SPLE2+D, 
are listed in Table 4. We do calculations for two differ- 
ent anisotropy radii, for = 00 (isotropic) and = i?cff 
(anisotropic). Since the observed stellar velocity disper- 

are completely negligible here. 



etc. 



13 



sion is so close to that already predicted from the lens 
models alone, we see that the final results barely differ 
from the results presented in §6.3, except that the error on 
Ho has decreased. In particular we find that the SPLE2 
models have nearly identical errors on Ho as the SPLEl+D 
models. Hence, both the lens and lens plus dynamics mod- 
els clearly prefer an isothermal mass distribution for Gl. 
In Table 5, we have listed the recovered point-image prop- 
erties from the best SPLEl-l-D model, assuming Vi = oo. 
Nearly identical results are found for = i?eff • The best- 
fit model of the Einstein Ring is shown in Fig. 5, along 
with the observed ring. Fig. 6 shows the critical and caus- 
tic curves and the time-delay surface. For comparison, the 
outer critical curve is also overplotted on the data on the 
left panel of Fig. 5. Note that the critical curve is consis- 
tent with going through the saddle points of the ring, as 
expected. This provides a consistency check of the model, 
independent of the assumption of ellipticity for the source. 

7. THE HUBBLE CONSTANT FROM B1608-I-656 

In this section we discuss the value of Ho from 
B1608-I-656 and its random and systematic errors in more 
detail, based on the class of models presented in §6. 

7.1. Models Using Only Lensing Constraints 

The SIE and SPLEl and SPLE2 models are mass mod- 
els based solely on the gravitational lensing constraints 
(i.e., no stellar dynamics) and give a value of Hq rang- 
ing from 71-74 km s~Hlpc"^ (Table 4). We adopt the 
value of Ho=74+;[° km s^^ Mpc^"^ as the best lensing-only 
determination, which includes the uncertainty on the den- 
sity slope of the Gl without any assumptions except for 
a prior on the slope of G2. We discuss the latter in more 
detail in §7.3. We note that including a prior on Gl sim- 
ilar to that on G2 does not change the result, since the 
slope of Gl is already found to be nearly isothermal, i.e., 
= 1.99t!];^!] (68% CL) without a prior. The lensing- 
only models therefore predict an isothermal density profile 
for Gl to within 10%. 

7.2. Models Using Lensing and Dynamical Constraints 

All other models we consider include the stellar velocity 
dispersion of Gl as a constraint. In §6, we showed that the 
simplest lens models that we consider, with two isother- 
mal galaxies, predict a stellar velocity dispersion close 
to the observed value (within l-tr for both the isotropic 
and strongly anisotropic models). This agreement shows 
that the slope of the density profile of Gl can not de- 
viate from isothermal too strongly. Indeed it is found 
that the slope of Gl is isothermal to within a few per- 
cent if we include the stellar velocity dispersion as an ad- 
ditional constraint. We finally adopt the median value 
of Ho=75lg km s~^ Mpc~^ as the best determination. A 
systematic error of ± 1 km s~^ Mpc~^ is estimated, due to 
the unknown anisotropy of Gl. We note that including the 
stellar dynamics has not significantly changed the value of 
Hq, but reduces its error considerably, to less than 10%. 
Similarly, the slope of Gl is = 2.03^5^:^4 ± 0-03 (68% 
CL), where we take the average value of the SPLEl+D 
models with = i?eff and = 00. The latter error is the 
uncertainty due to the unknown anisotropy of the stellar 
velocity elhpsoid of Gl (Table 4). 



7.3. Systematics 

Finally we discuss the dominant systematic errors that 
could potentially be left in the models. 

The results have not drastically changed since earlier 
models, but given the wealth of new constraints (e.g., the 
Einstein Ring, improved time delay measurements, stel- 
lar dynamics, etc.), the results have become more robust. 
The error on Ho has shrunk even though we have included 
more free parameters (e.g., external shear, varying slopes 
and galaxy positions, etc.) in the model. The smaller er- 
ror on Hq is due to improvements in the measurement of 
the time delays and the inclusion of the Einstein Ring and 
stellar dynamics. Support for the refined models is further 
given by the flux ratio of image D, which is recovered to 
within 22% (Tables 3 and 5), even though it has not been 
included as a constraint on the models. All previous mod- 
els had great diSiculty with this. Hence, at first glance, 
the results appear little affected by remaining systematics. 

Even so, we remain careful about assumptions, the most 
important of which is that for the prior on the density slope 
of G2 (§6.1). If we increase the width of the prior, we find 
that Ho increases somewhat (i.e., the density profile of G2 
becomes steeper), but also that its axial ratio decreases 
to very small values < 0.2 (PA remains similar to those 
in Table 4). Since G2 is most likely a disk galaxy (be- 
cause of the copious amounts of dust associated with it), 
a small axial ratio is not unlikely, but a value of ~0. 1-0.2 
appears unrealistic, especially since G2 has no discernible 
ellipticity in the F160W image (Fig. 1). If we take priors 
of qG2 =1.0±0.2 and loosen the prior on the density slope 
to 7q2=2.0±0.3 - both of which are also reasonable, given 
the observations - we find that the resulting value of Hq 
increases by at most a few percent. 

Next, we have assumed a spherical mass distribution 
for Gl in our stellar dynamical models. Since we measure 
a luminosity- weighted dispersion within a relatively large 
aperture (> Res), many assumptions about the phase- 
space density of stars have little effect on the final mea- 
sured stellar velocity dispersion. For example, anisotropy 
has little effect on Hq. A change of ±5% in the average 
dispersion would also affect Hq on a level of <3%. We 
therefore believe the assumption of spherical symmetry to 
be of little effect, given the current error on the observed 
stellar velocity dispersion of Gl (~14%). We also notice 
that the slope of Gl from lensing alone and from stellar dy- 
namics fully agree, giving no indication of large remaining 
systematic errors. 

As a check on the astrometry, dust correction and decon- 
volution procedure, we have used the Einstein Ring from 
the F814W observations as a constraint. We follow exactly 
the same procedure as for the F160W image. We find that 
the two lens galaxies lie slightly further apart than in the 
NICMOS images, as discussed previously. This suggests 
that a smaller value of Hq will be found. However, the 
best SPLEl-l-D model produces a somewhat steeper den- 
sity profile for Gl (by 2%), which compensates for the 
somewhat larger lens galaxy separation and leads to an 
identical value of Ho. Hence, it appears that inclusion of 
the ring "forces" the combined lens potential of Gl and G2 
to be slightly steeper than isothermal, leading to a robust 
value of Hq independent of which data set is used. 
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Fig. 6. — Left: Critical (thick) and caustic curves (thin) of the SPLEl+D model. The galaxy positions are indicated by stars, the images 
positions by open squares and the source by a closed square. Right: The contours indicate constant time delays starting at At=0 at image 
B and increasing in steps of 10 h^^ days. 



We therefore conclude that given the assumptions and 
priors - which are aU consistent with the observations of 
B1608+656 - remaining systematic errors on Hg probably 
add up to not more than about 5%, dominated by those 
from assumptions in the stellar dynamical model and the 
priors on G2. 

However, an important systematic error that could still 
affect the value of Hq is a mass-sheet degeneracy (Goren- 
stein, Shapiro, & Falco 1988) due to the small group 
around the lens system (Fassnacht et al. 2002a). Although 
our model reproduces the observed stellar velocity dis- 
persion well for an isothermal mass distribution, a strong 
mass-sheet could have led us to overestimate the enclosed 
galactic mass and therefore to underestimate the inferred 
density slope, given the observed velocity dispersion con- 
straint. The observed quasi-isothermality, however, does 
not necessarily argue against a mass-sheet since the mass 
density profile does not have to be exactly isothermal. We 
note also that the effect on the slope partially compen- 
sates the well known scaling of Ho with Kc for a given 
mass model. We defer a precise quantification of the con- 
vergence (kc) and shear of this group, and their effects on 
the lens system and Hq, to a forthcoming publication with 
additional data on the group (Fassnacht et al. in prep.), 
but believe it to be relatively small based on the "poor- 
ness" of the group and its small velocity dispersion. 

Finally, we note that the effects of the cosmological 
model are negligible. Changing Oa from 0.7 to 0.0, with 
fim = 0.3, increases Hq by 2%. A drastic change to 
(r^in = 1.0,17a = 0.0) decreases Hq by 7%. However, 
within the current limits on the cosmological model set 



by WMAP(Spergel et al. 2003), changes in Hq are < 1%. 

8. DISCUSSION & CONCLUSIONS 

We have presented significantly improved and refined 
mass models of the gravitational lens B1608-t-656 - com- 
pared with previous modeling efforts - with the aim of 
determining a value of the Hubble Constant, that is less 
affected by previously known systematics (e.g. radial mass 
profile, dust extinction, etc.). 

New constraints on the mass model include: (i) the stel- 
lar velocity dispersion of the dominant lens galaxy (Gl), 
as measured with ESI, (ii) the deconvolved Einstein Ring 
seen in the HST F160W and F814W images - the former 
of which is little affected by dust - corrected for the contri- 
bution from the lens galaxies, (iii) the extinction-corrected 
lens-galaxy centroids and structural parameters, the for- 
mer being one of the major uncertainties in previous lens 
models and (iv) recent improvements in the determination 
of the three independent time delays in this four-image lens 
system (Fassnacht et al. 2002b). 

Lens models have also been improved in allowing for 
many additional free parameters compared with previous 
modeling efforts (Koopmans & Fassnacht 1999; Fassnacht 
et al. 2002b), including the galaxy positions, the position 
angle of Gl, an external shear and the density slopes of Gl 
and G2. Some of these parameters are constrained with 
observational priors. The freedom in the lens model (up to 
22 free parameters) allows for a proper analysis of the er- 
ror on the inferred value of Hq, including all observational 
errors and correlations between free parameters. 

The improvements in the observations and lens models 
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have led to a Hubble Constant from B 1608+656 which is 
less affected by known systematic errors than previously. 
We obtain Ho=75t^ km s"! Mpc"^ (68% CL), with a ±5% 
error contributed by known systematic errors (for f2ni=0.3 
and Q/^=0.7). This value is higher than found from previ- 
ous models (e.g., Koopmans & Fassnacht 1999; Fassnacht 
et al. 2002b), predominantly because, after extinction cor- 
rection, the two lens galaxies are somewhat closer together 
than previously thought, their density profiles are slightly 
steeper than isothermal and external shear is included in 
the models. 

First, we find that the new lens models reproduce all 

radio-image constraints to within the errors. This agree- 
ment might imply that we have been too conservative in 
the error estimates and might also have slightly overesti- 
mated the final error on Hq. The flux ratio of image D 
is also reproduced to within 22%, a considerable improve- 
ment over previous models (e.g., Koopmans & Fassnacht 
1999; Fassnacht ot al. 2002b). Image D is the closest im- 
age to the dominant lens galaxy and might be affected by 
structure in the stellar (or dark-matter) mass distribution 
of lens galaxy Gl that is poorly represented by the global 
mass model. 

Second, the lens models reproduce the observed stellar 
velocity dispersion of Gl, and its mass profile is found to 
be nearly isothermal, = 2.031^14 ± 0.03 (68% CL). 
This result is not as straightforward as one might expect, 
since the inferred stellar velocity dispersion of Gl depends 
strongly on the mass that G2 contributes to the lens prop- 
erties. If the lens model predicted too large a mass for G2, 
the mass of Gl would decrease (because Icnsing tightly 
fixes the mass enclosed by the lensed images), as would the 
inferred stellar velocity dispersions. Similarly, if the mass 
is correct, but the density slope of Gl inferred from lensing 
alone is wrong, the inferred velocity dispersion would be 
incorrect as well. The fact that the lensing model (without 
dynamical constraints) already predicts a mass model for 
Gl that accurately reproduces the observed stellar velocity 
dispersions makes the models, in our opinion, very believ- 
able. The addition of the dynamical constraint therefore 
does not change the results from lensing alone, but signif- 
icantly tightens the error on the value Hq. 

Third, the lens models reproduce the boxy shape of the 
Einstein Ring both in the F160W and F814W bands well 
and give the same value of Hq irrespective of which of these 
datasets wc use as a constraint. In addition, the criti- 
cal curves (see Fig. 5) cross the Einstein Ring at the sad- 
dle points of its brightness distribution, as required (e.g., 
Blandford, Surpi, & Kundic 2001). 

We note that the resulting value Hq from B1608+656 
seems difficult to reconcile with Kochanek (2002, 2003), 
who finds that perfectly isothermal lens galaxies should 
lead on average to Hq = 48 ± 3 km s~^Mpc~^ (see also 
Kochanek & Schechter 2003). We find a significantly 
higher value of Hq even if the dominant lens galaxy (Gl) 
is assumed to be isothermal. One possible explanation is 
that the time-delay systems examined by Kochanek (2002, 
2003) have lens galaxies with luminous plus dark-matter 
mass profiles steeper than isothermal, enough to yield Hq 
in better agreement with that obtained from B1608-I-656. 
Although indeed - on average - isothermal models ap- 
pear to provide a good representation of the mass distribu- 



tion of early-type galaxies (e.g. Treu & Koopmans 2002a; 
Koopmans & Treu 2003, Rusin et al. 2003b), significant 
departures in several individual cases is not necessarily in- 
consistent with the observed spread of mass slopes of local 
early-type galaxies (e.g.. Fig. 1 in Gerhard et al. 2001; 
Bertin & Stiavelli 1993), especially at small radii where 
the stellar mass starts to dominate (e.g. Bertola et al. 
1993). The lens galaxy in PG1115+080 (Treu & Koop- 
mans 2002b) provides one such example of a lens-galaxy 
mass profile steeper than isothermal inside the Einstein 
radius. 

Conversely, a mass-sheet degeneracy could have led us 
to overestimate Hq from B1608+656, although the mass- 
sheet degeneracy is somewhat broken by the kinematic 
measurement. In this respect, we only note that a con- 
vergence of at least Kc~0.4 is needed to bring our value 
of Hq in agreement with 48 ± 3km Mpc~^, an unreal- 
istically high value, and wc defer further discussion to a 
forthcoming paper (Fassnacht et al. in prep.). 

Another viable alternative would be that the time-delay 
systems examined by Kochanek (2002, 2003) have not been 
modeled in enough detail yet and could also be affected by 
nearby clusters or groups as seen near some of these sys- 
tems (e.g., RXJ0911; Kneib, Cohen, & Hjorth 2000; Hjorth 
et al. 2002). A simple prediction is therefore the follow- 
ing: If these lens systems are isothermal and not affected 
by external perturbers, then the stellar velocity dispersion 
should be that predicted from a simple isothermal lens 
model (as for B1608-I-656). If, on the other hand, the ob- 
served stellar velocity dispersion is much higher that pre- 
dicted from isothermal models, this would be a clear sign 
of a steeper mass profile (as for PG1115-I-080). Hence, 
stellar kinematics of lens galaxies is a key measurement to 
further study this apparent inconsistency. 

Even so, the analysis of B1608-I-656 has demonstrated 
that when enough information is available, it is possible to 
reach accuracies of <10% (random) on the value of Hp (or 
<15% when including systematic errors) taking the radial 
mass profile properly into account. A similar uncertainty 
was obtained by applying a joint lensing and dynamical 
analysis to PGl 115+080 (Treu & Koopmans 2002b). This 
level of uncertainty is comparable to the accuracy achieved 
by the individual methods (SNae-Ia, Fundamental Plane, 
etc) that contributed to the final value of Hq from the HST 
Key-Project (Freedman et al. 2001). Hence, obtaining suf- 
ficiently good (lata including on the local environments 
- for a handful of gravitational- lens systems, seems to be a 
promising way to break the 10% limit on the uncertainty 
on the value of Hq . 

Finally, the analyses of B1608+656 and PGl 115+080 
have demonstrated that there are no simple ways to reach 
this goal and it does not suffice to assume a radial mass 
profile, without measuring it for each individual system. 
That approach probably does not lead to a reliable value 
of Hp. High-quality observations and detailed lens mass 
models - that allow for many degrees of freedom - are re- 
quired to obtain an accurate (and realistic) estimate of the 
Hubble Constant and its uncertainty. We are currently col- 
lecting spectroscopic data for additional lens systems with 
measured time-delays, with this aim in mind. In addition, 
optical and radio monitoring programs are continuing that 
should provide more (and/or improved) time delays in the 
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near future, and forthcoming observations with HST could 
provide improved constraints on the lens galaxies. 
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